home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / kcl / akcl / akcl1615.lha / mp / mp2.c < prev    next >
C/C++ Source or Header  |  1991-02-13  |  10KB  |  620 lines
  1.  
  2. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  3. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  4. /*~                                    ~*/
  5. /*~               OPERATIONS DE BASE (NOYAU)            ~*/
  6. /*~             Functions which can be efficient in plain C             ~*/
  7. /*~                                    ~*/
  8. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  9. /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
  10.  
  11.  
  12. #include "config.h"
  13. #include "genpari.h"
  14. #include "arith.h"
  15.  
  16. /* -2147483648 */
  17.  
  18. ulong MOST_NEGS[3]={0x01ff0003, 0xff000003,1<<31};
  19.  
  20. /* +2147483648 */
  21.  
  22. ulong ABS_MOST_NEGS[3]={0x01ff0003, 0x01000003,1<<31};
  23.  
  24.  
  25. GEN stoi(x)
  26.      long x;
  27. {
  28.   GEN y;
  29.   
  30.   if(!x) return gzero;
  31.   y=cgeti(3);
  32.   if(x>0) {y[1]=0x1000003;y[2]=x;}
  33.   else{y[1]=0xff000003;y[2]= -x;}
  34.   return y;
  35. }
  36.  
  37.  
  38. GEN cgetg(x,y)
  39.      long x,y;
  40. {
  41.   unsigned long p1;
  42.   GEN z;
  43.   
  44.   p1=avma-(((unsigned short)x)<<2);if(p1<bot) err(errpile);
  45.   avma=p1;z=(GEN)p1;z[0]=0x10000+x+(y<<24);
  46.   return z;
  47. }
  48.  
  49.  
  50. GEN cgeti(x)
  51.      long x;
  52. {
  53.   unsigned long p1;
  54.   GEN z;
  55.   
  56.   p1=avma-4*x;if(p1<bot) err(errpile);
  57.   avma=p1;z=(GEN)p1;z[0]=0x1010000+x;
  58.   return z;
  59. }
  60.  
  61.  
  62. GEN icopy(x)
  63.      GEN x;
  64. {
  65.   GEN y;
  66.   long lx=lgef(x),i;
  67.   
  68.   y=cgeti(lx);
  69.   for(i=1;i<lx;i++) y[i]=x[i];
  70.   return y;
  71. }
  72.  
  73.  
  74. GEN negi(x)
  75.      GEN x;
  76. {
  77.   long s=signe(x);
  78.   GEN y;
  79.   
  80.   if(!s) return gzero;
  81.   y=icopy(x);setsigne(y,-s);
  82.   return y;
  83. }
  84.  
  85.  
  86. GEN absi(x)
  87.      GEN x;
  88. {
  89.   GEN y;
  90.   long s=signe(x);
  91.   
  92.   if(!s) return gzero;
  93.   y=icopy(x);setsigne(y,1);return y;
  94. }
  95.  
  96.  
  97. long itos(x)
  98.      GEN x;
  99. {
  100.   long s=signe(x),p2;
  101.   unsigned long p1;
  102.   
  103.   if(!s) return 0;
  104.   if(lgef(x)>3) err(affer2);
  105.   p1=x[2];if(p1>=0x80000000) err(affer2);
  106.   p2=(s>0)?p1:(-((long)p1));return p2;
  107. }
  108.  
  109.  
  110. void affsi(s,x)
  111.      long s;
  112.      GEN x;
  113. {
  114.   long lx;
  115.   
  116.   if(!s) {x[1]=2;return;}
  117.   lx=lg(x);if(lx<3) err(affer1);
  118.   if(s>0) {x[1]=0x1000003;x[2]=s;}
  119.   else { s = -s;
  120.      if (s < 0) /* s = -2^31 */
  121.        { if(lx<4) err(affer1);
  122.          x[1]=0xff000004;
  123.          x[2]= 0;
  124.          x[3]= 1;
  125.        }
  126.        else 
  127.          {x[1]=0xff000003;x[2]= s;}
  128.   }
  129. }
  130.  
  131.  
  132. void affii(x,y)
  133.      GEN x,y;
  134. {
  135.   long lx=lgef(x),i;
  136.   
  137.   if(x==y) return;
  138.   if(lg(y)<lx) err(affer3);
  139.   for(i=1;i<lx;i++) y[i]=x[i];
  140. }
  141.  
  142.  
  143. GEN shifts(x,y)
  144.      long x,y;
  145. {
  146.   long t[3];
  147.   
  148.   if(!x) return gzero;
  149.   t[0]=0x1010003;
  150.   if(x>0) {t[1]=0x1000003;t[2]=x;}
  151.   else {t[1]=0xff000003;t[2]= -x;}
  152.   return shifti(t,y);
  153. }
  154.  
  155.  
  156. GEN shifti(x,n)
  157.      GEN x;
  158.      long n;
  159. { 
  160.   long lx=lgef(x),i,s=signe(x),d,m,p1,p2,k;
  161.   GEN y; TEMPVARS
  162.   ulong hiremainder;
  163.   
  164.   if(!s) return gzero;
  165.   if(!n) return icopy(x);
  166.   if(n>0)
  167.     {
  168.       d=n>>5;m=n&31;
  169.       if(m)
  170.     {
  171.       p1=shiftl(x[2],m);p2=hiremainder;k=0;
  172.       if(p2)
  173.         {
  174.           y=cgeti(lx+d+1);for(i=lx+1;i<=lx+d;i++) y[i]=0;
  175.           for(i=lx;i>=4;i--) {y[i]=shiftl(x[i-1],m)+k;k=hiremainder;}
  176.           y[3]=p1+k;y[2]=p2;
  177.         }
  178.       else
  179.         {
  180.           y=cgeti(lx+d);for(i=lx;i<lx+d;i++) y[i]=0;
  181.           for(i=lx-1;i>=3;i--) {y[i]=shiftl(x[i],m)+k;k=hiremainder;}
  182.           y[2]=p1+k;
  183.         }
  184.     }
  185.       else
  186.     {
  187.       y=cgeti(lx+d);for(i=lx;i<lx+d;i++) y[i]=0;
  188.       for(i=lx-1;i>=2;i--) y[i]=x[i];
  189.     }
  190.     }
  191.   else
  192.     {
  193.       n= -n;d=n>>5;m=n&31;if(lx<d+3) return gzero;
  194.       if(!m)
  195.     {
  196.       y=cgeti(lx-d);for(i=2;i<lx-d;i++) y[i]=x[i];
  197.     }
  198.       else 
  199.     {
  200.       m=32-m;d++;
  201.       p1=shiftl(x[2],m);
  202.       if(hiremainder)
  203.         {
  204.           y=cgeti(lx-d+1);y[2]=hiremainder;
  205.           for(i=3;i<=lx-d;i++)
  206.         {
  207.           p2=shiftl(x[i],m);y[i]=p1+hiremainder;p1=p2;
  208.         }
  209.         }
  210.       else
  211.         {
  212.           if(lx==d+2) return gzero;
  213.           y=cgeti(lx-d);
  214.           for(i=3;i<=lx-d;i++) 
  215.         {
  216.           p2=shiftl(x[i],m);y[i-1]=p1+hiremainder;p1=p2;
  217.         }
  218.         }
  219.     }
  220.     }   
  221.   y[1]=y[0];setsigne(y,s);return y;
  222. }
  223.  
  224.  
  225. GEN mptrunc(x)
  226.      GEN x;
  227. {
  228.   long e,i,s,lx=lg(x),p1,p2,m;
  229.   unsigned long d; TEMPVARS
  230.   GEN y;ulong hiremainder;
  231.   
  232.   if(typ(x)==1) return icopy(x);
  233.   s=signe(x);if(!s) return gzero;
  234.   e=expo(x);if(e<0) return gzero;
  235.   d=e>>5;m=e&31;if(d>=lx-2) err(truer2);
  236.   y=cgeti(d+3);y[1]=y[0];setsigne(y,s);
  237.   if(m==31) for(i=2;i<=d+2;i++) y[i]=x[i];
  238.   else
  239.     {
  240.       m++;p1=0;
  241.       for(i=2;i<=d+2;i++)
  242.     {
  243.       p2=shiftl(x[i],m);y[i]=hiremainder+p1;p1=p2;
  244.     }
  245.     }
  246.   return y;
  247. }
  248.  
  249.  
  250. GEN mpent(x)
  251.      GEN x;
  252. {
  253.   long e,i,lx=lg(x),m,f,p1,p2;
  254.   unsigned long d;ulong hiremainder;
  255.   GEN y,z; TEMPVARS
  256.   
  257.   if(typ(x)==1) return icopy(x);
  258.   if(signe(x)>=0) return mptrunc(x);
  259.   e=expo(x);if(e<0) {y=cgeti(3);y[2]=1;y[1]=0xff000003;return y;}
  260.   d=e>>5;m=e&31;if(d>=lx-2) err(truer2);
  261.   y=cgeti(d+3);y[1]=0xff000003+d;
  262.   if(m==31) 
  263.     {
  264.       for(i=2;i<=d+2;i++) y[i]=x[i];
  265.       while((i<lx)&&(!x[i])) i++;
  266.       f=(i<lx);
  267.     }    
  268.   else
  269.     {
  270.       m++;p1=0;
  271.       for(i=2;i<=d+2;i++)
  272.     {
  273.       p2=shiftl(x[i],m);y[i]=hiremainder+p1;p1=p2;
  274.     }
  275.       if(p1) f=1;
  276.       else
  277.     {
  278.       while((i<lx)&&(!x[i])) i++;
  279.       f=(i<lx);
  280.     }
  281.     }
  282.   if(f)
  283.     {
  284.       for(i=d+2;(i>=2)&&(y[i]==0xffffffff);i--) y[i]=0;
  285.       if(i>=2) y[i]++;
  286.       else
  287.     {
  288.       z=y;y=cgeti(1);*y=(*z)+1;y[1]=z[1]+1;
  289.     }
  290.     }
  291.   return y;
  292. }
  293.  
  294.  
  295. int cmpsi(x,y)
  296.      long x;
  297.      GEN y;
  298. {
  299.   ulong p;
  300.   
  301.   if(!x) return -signe(y);
  302.   if(x>0)
  303.     {
  304.       if(signe(y)<=0) return 1;
  305.       if(lgef(y)>3) return -1;
  306.       p=y[2];if(p==x) return 0;
  307.       return (p<(ulong)x) ? 1 : -1;
  308.     }
  309.   else
  310.     {  /* x <= 0 */
  311.       if(signe(y)>=0) return -1;
  312.       if(lgef(y)>3)
  313.     { if (-x < 0)
  314.         { /* x = -2^31 */
  315.           if (lgef(y)==4 &&
  316.           y[2] == 0 &&
  317.           y[3] == 1)
  318.         return 0;
  319.         else
  320.           return 1;}}
  321.       p=y[2];if(p== -x) return 0;
  322.       return (p<(ulong)(-x)) ? -1 : 1;
  323.     }
  324. }
  325.  
  326.  
  327. int cmpii(x,y)
  328.      GEN x,y;
  329. {
  330.   long sx=signe(x),sy=signe(y),lx,ly,i;
  331.   
  332.   if(sx<sy) return -1;
  333.   if(sx>sy) return 1;
  334.   if(!sx) return 0;
  335.   lx=lgef(x);ly=lgef(y);
  336.   if(lx>ly) return sx;
  337.   if(lx<ly) return -sx;
  338.   for(i=2;(i<lx)&&(x[i]==y[i]);i++);
  339.   if(i==lx) return 0;
  340.   return ((ulong)x[i]>(ulong)y[i]) ? sx : -sx;
  341. }
  342.  
  343.  
  344. GEN addss(x,y)
  345.      long x,y;
  346. {
  347.   long t[3];
  348.   
  349.   if(!x) return stoi(y);
  350.   t[0]=0x1010003;
  351.   if(x>0) {t[1]=0x1000003;t[2]=x;} else {t[1]=0xff000003;t[2]= -x;}
  352.   return addsi(y,t);
  353. }
  354.  
  355.  
  356. GEN subii(x,y)
  357.      GEN x,y;
  358. {
  359.   long s=signe(y);
  360.   GEN z;
  361.   
  362.   if(x==y) return gzero;
  363.   setsigne(y,-s);z=addii(x,y);setsigne(y,s);
  364.   return z;
  365. }
  366.  
  367.  
  368. GEN subsi(x,y)
  369.      long x;
  370.      GEN y;
  371. {
  372.   long s=signe(y);
  373.   GEN z;
  374.   
  375.   setsigne(y,-s);z=addsi(x,y);setsigne(y,s);return z;
  376. }
  377.  
  378.  
  379. GEN subss(x,y)
  380.      long x,y;
  381. {
  382.   if (y == (1<<31))
  383.     return addsi(x,ABS_MOST_NEGS);
  384.   return addss(-y,x);
  385. }
  386.  
  387.  
  388. GEN convi(x)
  389.      GEN x;
  390. {
  391.   long lx,av=avma,lz;
  392.   GEN z,p1,p2;  
  393.   
  394.   if(!signe(x))
  395.     {
  396.       z=cgeti(3);z[1]= -1;z[2]=0;avma=av;return z+3;
  397.     }
  398.   p1=absi(x);lx=lgef(p1);lz=((lx-2)*15)/14+3;z=cgeti(lz);z[1]= -1;
  399.   for(p2=z+2;signe(p1);p2++) *p2=divisii(p1,1000000000,p1);
  400.   avma=av;return p2;
  401. }
  402.  
  403.  
  404.  
  405. void mulsii(x,y,z)
  406.      long x;
  407.      GEN y,z;
  408. {
  409.   long av=avma;
  410.   GEN p1;
  411.   
  412.   p1=mulsi(x,y);affii(p1,z);avma=av;
  413. }
  414.  
  415.  
  416. void addsii(x,y,z)
  417.      long x;
  418.      GEN y,z;
  419. {
  420.   long av=avma;
  421.   GEN p1;
  422.   
  423.   p1=addsi(x,y);affii(p1,z);avma=av;
  424. }
  425.  
  426.  
  427. long divisii(x,y,z)
  428.      long y;
  429.      GEN x,z;
  430. {
  431.   long av=avma,k;
  432.   GEN p1;
  433.   
  434.   p1=divis(x,y);affii(p1,z);avma=av;
  435.   k=hiremainder;return k;
  436. }
  437.  
  438.  
  439. long vals(x)
  440.      long x;
  441. {
  442.   unsigned short int y,z;
  443.   int s;
  444.  
  445.   if(!x) return -1;
  446.   y=x;if(!y) {s=16;y=((ulong)x)>>16;} else s=0;
  447.   z=y&255;if(!z) {s+=8;z=y>>8;}
  448.   y=z&15;if(!y) {s+=4;y=z>>4;}
  449.   z=y&3;if(!z) {s+=2;z=y>>2;}
  450.   return (z&1) ? s : s+1;
  451. }
  452.  
  453.  
  454. long vali(x)
  455.      GEN x;
  456. {
  457.   long i,lx=lgef(x);
  458.   
  459.   if(!signe(x)) return -1;
  460.   for(i=lx-1;(i>=2)&&(!x[i]);i--);
  461.   return ((lx-1-i)<<5)+vals(x[i]);
  462. }
  463.  
  464. GEN dvmdss(x,y,z)
  465.      long x,y;
  466.      GEN *z;
  467. {
  468.   GEN p1;
  469.  
  470.   p1=divss(x,y);*z=stoi(hiremainder);
  471.   return p1;
  472. }
  473.  
  474.  
  475. GEN dvmdsi(x,y,z)
  476.      long x;
  477.      GEN y,*z;
  478. {
  479.   GEN p1;
  480.   p1=divsi(x,y);*z=stoi(hiremainder);
  481.   return p1;
  482. }
  483.  
  484.  
  485. GEN dvmdis(x,y,z)
  486.      long y;
  487.      GEN x,*z;
  488. {
  489.   GEN p1;
  490.   p1=divis(x,y);*z=stoi(hiremainder);
  491.   return p1;
  492. }
  493.  
  494.  
  495. GEN ressi(x,y)
  496.      long x;
  497.      GEN y;
  498. { 
  499.   divsi(x,y);return stoi(hiremainder);
  500. }
  501.  
  502.  
  503. GEN modsi(x,y)
  504.      long x;
  505.      GEN y;
  506. {
  507.   long s;
  508.   GEN p1;
  509.   
  510.   divsi(x,y);
  511.   if(!hiremainder) return gzero;
  512.   if(x>0) return stoi(hiremainder);
  513.   else
  514.     {
  515.       s=signe(y);setsigne(y,1);p1=addsi(hiremainder,y);
  516.       setsigne(y,s);return p1;
  517.     }
  518. }
  519.  
  520.  
  521. GEN modis(x,y)
  522.      long y;
  523.      GEN x;
  524. { 
  525.   divis(x,y);if(!hiremainder) return gzero;
  526.   return (signe(x)>0) ? stoi(hiremainder) : stoi(abs(y)+hiremainder);
  527. }
  528.  
  529.  
  530. GEN resis(x,y)
  531.      long y;
  532.      GEN x;
  533. { 
  534.   divis(x,y);return stoi(hiremainder);
  535. }
  536.  
  537.  
  538. GEN modii(x,y)
  539.      GEN x,y;
  540. {
  541.   long av=avma,tetpil;
  542.   GEN p1;
  543.  
  544.   p1=dvmdii(x,y,-1);
  545.   if(signe(p1)>=0) return p1;
  546.   tetpil=avma;p1=(signe(y)>0) ? addii(p1,y) : subii(p1,y);
  547.   return gerepile(av,tetpil,p1);
  548. }
  549.  
  550.  
  551. mpdivis(x,y,z)
  552.      GEN x,y,z;
  553. {
  554.   long av=avma;
  555.   GEN p1,p2;
  556.  
  557.   p1=dvmdii(x,y,&p2);
  558.   if(signe(p2)) {avma=av;return 0;}
  559.   affii(p1,z);avma=av;return 1;
  560. }
  561.  
  562.  
  563. divise(x,y)
  564.      GEN x,y;
  565. {
  566.   long av=avma;
  567.   GEN p1;
  568.  
  569.   p1=dvmdii(x,y,-1);avma=av;
  570.   return signe(p1) ? 0 : 1;
  571. }
  572.  
  573.  
  574. GEN gerepile(l,p,q)
  575.     GEN l,p,q;
  576.  
  577. {
  578.   long av,declg,tl;
  579.   GEN ll,pp,l1,l2,l3;
  580.  
  581.   declg=(long)l-(long)p;if(declg<=0) return q;
  582.   for(ll=l,pp=p;pp>(GEN)avma;) *--ll= *--pp;
  583.   av=(long)ll;
  584.   while((ll<l)||((ll==l)&&(long)q))
  585.   {
  586.     l2=ll+lontyp[tl=typ(ll)];
  587.     if(tl==10) {l3=ll+lgef(ll);ll+=lg(ll);if(l3>ll) l3=l2;}
  588.     else {ll+=lg(ll);l3=ll;} 
  589.     for(;l2<l3;l2++) 
  590. /*    for(;l2<ll;l2++) */
  591.       {
  592.     l1=(GEN)(*l2);
  593.     if((l1<l)&&(l1>=(GEN)avma))
  594.       {
  595.         if(l1<p) *l2+=declg;
  596.         else
  597.           if(ll<=l) err(gerper);
  598.       }
  599.       }
  600.   }
  601.   if((!((long)q))||((q<p)&&(q>=(GEN)avma)))
  602.   {
  603.     avma=av;return q+(declg>>2);
  604.   }
  605.   else {avma=av;return q;}
  606. }
  607.  
  608.  
  609. void cgiv(x)
  610.      GEN x;
  611. {
  612.   long p;
  613.  
  614.   if((p=pere(x))==255) return;
  615.   if((x!=(GEN)avma)||(p>1)) {setpere(x,p-1);return;}
  616.   do x+=lg(x);while(!pere(x));
  617.   avma=(long)x;
  618.   return;
  619. }
  620.